# 01_qc_preprocessing.R
# Seurat单细胞数据QC和预处理

# 输入: 原始单细胞数据
# 输出: seurat_obj_qc.rds, QC指标可视化

library(Seurat)
library(ggplot2)
library(dplyr)
library(patchwork)

set.seed(123)

cat("=== Seurat单细胞数据QC和预处理开始 ===\n")

# 创建输出目录
if (!dir.exists("output")) {
  dir.create("output", recursive = TRUE)
}
if (!dir.exists("figures")) {
  dir.create("figures", recursive = TRUE)
}

# 检查数据文件是否存在
data_dir <- "data"
if (!dir.exists(data_dir)) {
  stop("错误: 找不到data目录，请确保单细胞数据文件在data/目录下")
}

# 读取数据
cat("读取单细胞数据...\n")

# 方法1: 如果是从10X Genomics Cell Ranger输出读取
if (file.exists(file.path(data_dir, "matrix.mtx")) && 
    file.exists(file.path(data_dir, "features.tsv")) && 
    file.exists(file.path(data_dir, "barcodes.tsv"))) {
  
  cat("检测到10X Genomics格式数据，正在读取...\n")
  pbmc.data <- Read10X(data.dir = data_dir)
  seurat_obj <- CreateSeuratObject(counts = pbmc.data, 
                                  project = "pbmc3k", 
                                  min.cells = 3, 
                                  min.features = 200)
  
} else if (file.exists(file.path(data_dir, "pbmc3k_filtered_gene_bc_matrices.tar.gz"))) {
  
  # 方法2: 如果是从教程的压缩文件读取
  cat("检测到PBMC3K教程数据，正在解压和读取...\n")
  untar(file.path(data_dir, "pbmc3k_filtered_gene_bc_matrices.tar.gz"), exdir = data_dir)
  pbmc.data <- Read10X(data.dir = file.path(data_dir, "filtered_gene_bc_matrices/hg19"))
  seurat_obj <- CreateSeuratObject(counts = pbmc.data, 
                                  project = "pbmc3k", 
                                  min.cells = 3, 
                                  min.features = 200)
  
} else {
  
  # 方法3: 尝试直接加载内置数据或示例数据
  cat("尝试加载示例数据...\n")
  # 这里可以使用内置的pbmc3k数据
  if (!file.exists("output/pbmc3k_raw.rds")) {
    # 下载或使用示例数据
    pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
    seurat_obj <- CreateSeuratObject(counts = pbmc.data, 
                                    project = "pbmc3k", 
                                    min.cells = 3, 
                                    min.features = 200)
    saveRDS(seurat_obj, "output/pbmc3k_raw.rds")
  } else {
    seurat_obj <- readRDS("output/pbmc3k_raw.rds")
  }
}

cat("原始数据维度:", dim(seurat_obj), "\n")
cat("细胞数:", ncol(seurat_obj), "\n")
cat("基因数:", nrow(seurat_obj), "\n")

# 计算线粒体基因百分比
cat("计算QC指标...\n")
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

# 计算核糖体基因百分比（可选）
seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# 计算红细胞基因百分比（可选，针对PBMC）
seurat_obj[["percent.hb"]] <- PercentageFeatureSet(seurat_obj, pattern = "^HB[^(P)]")

# 查看QC指标摘要
cat("\nQC指标摘要:\n")
print(summary(seurat_obj@meta.data$nFeature_RNA))
print(summary(seurat_obj@meta.data$nCount_RNA))
print(summary(seurat_obj@meta.data$percent.mt))

# 可视化QC指标
cat("生成QC指标可视化...\n")

# 创建QC指标的数据框用于绘图
qc_data <- seurat_obj@meta.data

# 1. 特征数（基因数）分布
p1 <- ggplot(qc_data, aes(x = nFeature_RNA)) +
  geom_histogram(bins = 50, fill = "lightblue", color = "black") +
  geom_vline(xintercept = c(200, 2500), linetype = "dashed", color = "red") +
  labs(title = "Distribution of Genes per Cell",
       x = "Number of Genes",
       y = "Number of Cells") +
  theme_minimal()

# 2. UMI计数分布
p2 <- ggplot(qc_data, aes(x = nCount_RNA)) +
  geom_histogram(bins = 50, fill = "lightgreen", color = "black") +
  geom_vline(xintercept = c(500, 10000), linetype = "dashed", color = "red") +
  labs(title = "Distribution of UMIs per Cell",
       x = "Number of UMIs",
       y = "Number of Cells") +
  theme_minimal()

# 3. 线粒体基因百分比分布
p3 <- ggplot(qc_data, aes(x = percent.mt)) +
  geom_histogram(bins = 50, fill = "lightcoral", color = "black") +
  geom_vline(xintercept = 5, linetype = "dashed", color = "red") +
  labs(title = "Distribution of Mitochondrial Gene Percentage",
       x = "Mitochondrial Gene Percentage (%)",
       y = "Number of Cells") +
  theme_minimal()

# 组合QC图
qc_histograms <- p1 / p2 / p3
ggsave("output/qc_histograms.png", qc_histograms, width = 8, height = 10, dpi = 300)

# 散点图展示QC指标之间的关系
p4 <- ggplot(qc_data, aes(x = nCount_RNA, y = nFeature_RNA, color = percent.mt)) +
  geom_point(alpha = 0.7, size = 0.5) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "UMIs vs Genes colored by MT%",
       x = "Number of UMIs",
       y = "Number of Genes",
       color = "MT%") +
  theme_minimal()

p5 <- ggplot(qc_data, aes(x = nCount_RNA, y = percent.mt)) +
  geom_point(alpha = 0.7, size = 0.5, color = "darkblue") +
  geom_hline(yintercept = 5, linetype = "dashed", color = "red") +
  labs(title = "UMIs vs Mitochondrial Percentage",
       x = "Number of UMIs",
       y = "Mitochondrial Gene Percentage (%)") +
  theme_minimal()

p6 <- ggplot(qc_data, aes(x = nFeature_RNA, y = percent.mt)) +
  geom_point(alpha = 0.7, size = 0.5, color = "darkgreen") +
  geom_hline(yintercept = 5, linetype = "dashed", color = "red") +
  labs(title = "Genes vs Mitochondrial Percentage",
       x = "Number of Genes",
       y = "Mitochondrial Gene Percentage (%)") +
  theme_minimal()

qc_scatter <- (p4 | p5 | p6)
ggsave("output/qc_scatter.png", qc_scatter, width = 15, height = 5, dpi = 300)

# 使用Seurat内置函数可视化QC指标
qc_vln <- VlnPlot(seurat_obj, 
                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
                  ncol = 3)
ggsave("output/qc_violin.png", qc_vln, width = 12, height = 5, dpi = 300)

# 应用QC过滤
cat("应用QC过滤...\n")

# 保存过滤前的细胞数
cells_before <- ncol(seurat_obj)

# 设置QC阈值（可根据数据调整）
min_genes <- 200      # 最少基因数
max_genes <- 2500     # 最多基因数
min_umis <- 500       # 最少UMI数
max_umis <- 10000     # 最多UMI数
max_mt_percent <- 5   # 最大线粒体基因百分比

# 应用过滤
seurat_obj <- subset(seurat_obj, 
                     subset = nFeature_RNA > min_genes & 
                              nFeature_RNA < max_genes & 
                              nCount_RNA > min_umis & 
                              nCount_RNA < max_umis & 
                              percent.mt < max_mt_percent)

# 保存过滤后的细胞数
cells_after <- ncol(seurat_obj)

cat("过滤前细胞数:", cells_before, "\n")
cat("过滤后细胞数:", cells_after, "\n")
cat("过滤掉的细胞数:", cells_before - cells_after, "\n")
cat("保留细胞百分比:", round(cells_after/cells_before * 100, 2), "%\n")

# 可视化过滤后的QC指标
qc_vln_filtered <- VlnPlot(seurat_obj, 
                          features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
                          ncol = 3)
ggsave("output/qc_violin_filtered.png", qc_vln_filtered, width = 12, height = 5, dpi = 300)

# 保存过滤后的对象
cat("保存QC后的Seurat对象...\n")
saveRDS(seurat_obj, "output/seurat_obj_qc.rds")

# 生成QC报告
qc_report <- data.frame(
  Metric = c("Initial_Cells", "After_QC_Cells", "Removed_Cells", "Retention_Rate",
             "Min_Genes", "Max_Genes", "Min_UMIs", "Max_UMIs", "Max_MT_Percent"),
  Value = c(cells_before, cells_after, cells_before - cells_after, 
           round(cells_after/cells_before * 100, 2),
           min_genes, max_genes, min_umis, max_umis, max_mt_percent)
)

write.csv(qc_report, "output/qc_report.csv", row.names = FALSE)

# 生成过滤前后的对比图
filter_comparison <- data.frame(
  Stage = rep(c("Before QC", "After QC"), each = 3),
  Metric = rep(c("Genes", "UMIs", "MT%"), 2),
  Value = c(median(seurat_obj@meta.data$nFeature_RNA), 
           median(seurat_obj@meta.data$nCount_RNA), 
           median(seurat_obj@meta.data$percent.mt),
           median(qc_data$nFeature_RNA), 
           median(qc_data$nCount_RNA), 
           median(qc_data$percent.mt))
)

p_comparison <- ggplot(filter_comparison, aes(x = Stage, y = Value, fill = Metric)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~Metric, scales = "free_y") +
  labs(title = "QC Metrics Before and After Filtering",
       y = "Median Value") +
  theme_minimal() +
  scale_fill_brewer(palette = "Set1")

ggsave("output/qc_comparison.png", p_comparison, width = 10, height = 6, dpi = 300)

cat("=== QC和预处理完成 ===\n")
cat("输出文件:\n")
cat("- output/seurat_obj_qc.rds: QC后的Seurat对象\n")
cat("- output/qc_histograms.png: QC指标分布直方图\n")
cat("- output/qc_scatter.png: QC指标散点图\n")
cat("- output/qc_violin*.png: QC指标小提琴图\n")
cat("- output/qc_comparison.png: 过滤前后对比图\n")
cat("- output/qc_report.csv: QC报告\n")

cat("\nQC过滤统计:\n")
cat("- 初始细胞数:", cells_before, "\n")
cat("- 过滤后细胞数:", cells_after, "\n")
cat("- 细胞保留率:", round(cells_after/cells_before * 100, 2), "%\n")
cat("- 过滤阈值: nFeature_RNA ∈ [", min_genes, ", ", max_genes, "], ", 
    "nCount_RNA ∈ [", min_umis, ", ", max_umis, "], ", 
    "percent.mt < ", max_mt_percent, "\n", sep = "")